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Abstract: Splines are a popular and attractive way of smoothing noisy data. Computing splines in- 
volves minimizing a functional which is a linear combination of a fitting term and a regularization term. 
The former is classically computed using a (weighted) L2 norm while the latter ensures smoothness. 
Thus, when dealing with grid data, the optimization can be solved very efficiently using the DCT. In 
this work we propose to replace the L2 norm in the fitting term with an LI norm, leading to automatic 
robustness to outliers. To solve the resulting minimization problem we propose an extremely simple 
and efficient numerical scheme based on split-Bregman iteration combined with DCT. Experimental 
validation shows the high-quality results obtained in short processing times. 
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1. Introduction 

Smoothing a dataset consists in finding an approximating function that captures important patterns in the 
data, while disregarding noise or other fine-scale structures. Let y G W llX "' xnrn — » C be an m-dimensional 
discrete signal, where rij (1 < j < d) is the domain of y along the j-th dimension. We can model y by 

y = y + r, (1) 

where r represents some noise and y is a smooth function. A very common regularization choice is to enforce 
C 2 continuity, in which case y is called a (cubic) spline. Smoothing y relies upon finding the best estimate 
of y under the proper smoothness and noise assumptions. We can approximate y by minimizing an objective 
functional 

F(z) = R y (z) + sP(z), (2) 

where R y (z) is a data fitting term, defined by the distribution of r, P(z) is a regularization term, and s is a 
scalar that determines the balance between both terms. Such scalar parameter can be automatically derived 
using Bayesian or MDL techniques, as will be later shown. 

For clarity, we will describe in depth the case m = 1, where we have n = n\ samples. We extend the 
results later to the general m-dimensional case. 

Let us begin by explaining the smoothing term. The C 2 continuity requirement leads to define P(^) = 
||Dz||2, D being a discrete second-order differential operator, defined Vi, 2 < i < n — 1, by 

D hi~l = h i -i{h i - 1 +h i )i ( 3 ) 

= *A> ( 4 ) 

= hiihi-^+hi)! ( 5 ) 

where hi represents the step, or sampling rate, between yi and Assuming repeating border elements, 

that is, 2/0 = Vi and y n+1 = y n , gives D hl = -D ia = -l/h\, and D n>n -i = -Ai,n = -l/ ft n- 



l- 



Regarding the fitting term, the classical assumption is that the noise r in Equation (1) has Gaussian 

\l 



distribution with zero mean and unknown variance, which leads to setting R y (z) = \\z — y\\^. Smoothing 



then can be formulated as the least-squares regression 



y = argmin \\z - yf 2 + s \\Dzf 2 , 



*This work was partially done while the authors were with the Department of Electrical and Computer Engineering, Uni- 
versity of Minnesota. 
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For clarity, we call the estimate y obtained with this method an L2 spline. 

It is a well known fact that least squares estimates for regression models are highly non-robust to outliers. 
Although there is no agreement on a universal and formal definition of an outlier, it is usually regarded 
as an observation that does not follow the patterns in the data. Notice that smoothing should produce an 
estimate y taking into account only important patterns, that is, the inliers, in the data. In this sense, the L2 
formulation cannot correctly handle outliers by itself. 

In order to solve this problem, we propose to take a different assumption on the distribution of the noise 
r in Equation (1). By choosing a distribution with fatter tails than the Gaussian distribution, the derived 
estimator will correctly handle outliers. We thus assume that r follows a Laplace distribution with zero mean 
and unknown scale parameter, a common practice in other problems as we will further discuss below. This 
leads to the fitting term K y (z) = \\z — y^ and the regression then becomes 

y = argmin \\z - y\\ x + s \\Dz\\ 2 2 . 

z 

We call the estimate y obtained with this formulation an LI spline. 

Let us point out that the use of LI fitting terms for solving inverse problems is not new. For example, 
in 2001, Nikolova proved the theoretic pertinence of using LI fitting terms for image denoising [11]. Other 
interesting works have addressed this approach for total variation image denoising [1,2,5,12] or total variation 
optical flow [14,16,21] (this robustness-to-outliers type of ideas was previously introduced in the context of 
optical flow by Black and Anandan [3]). Recall that the total- variation regularization term involves first-order 
derivatives, while the proposed LI splines, on the other hand, involve second-order derivatives. 

Regularly sampled signals are extremely common in practice, and their analysis becomes easier and faster. 
In particular, we follow the most common choice when dealing with discrete m-dimensional data, which is 
assuming a "rectangular" Cartesian sampling pattern. When the sampling is isotropic, i.e., "square," we 
refer to this type of data as grid data. 

We developed an iterative algorithm for computing LI splines, based on split-Bregman iteration [8], that 
is specially suited for the case of grid data. This algorithm is extremely fast, both in running time and 
in the number of iterations until convergence. It is also outstandingly simple, making the implementation 
completely straightforward. 

The remainder of the paper is structured as follows. In Section 2, we overview a fast algorithm to compute 
L2 splines and robust L2 splines, a modification of the least squares regression that allows to handle outliers. 
In Section 3, we present the proposed algorithm for computing LI splines. Then, in Section 4, we show results 
obtained with LI splines, systematically outperforming its L2 and robust L2 counterparts in the presence of 
outliers. We also show that the proposed computational algorithm is very efficient. Finally, in Section 5 we 
provide some concluding remarks. 



2. Smoothing splines 

As aforementioned, the classical assumption is that the noise r in Equation (1) follows a Gaussian distribution 
with zero mean and unknown variance. This leads to solve the least-squares regression 



y = argmin \\z - y\\\ + s \\Dzf 2 , 



(6) 



Since both terms are different iable, we obtain 

y=(I + sD T D)- 1 y. 



(7) 



Garcia proposed a very efficient method for dealing with regularly sampled data [7]. Assuming that the 
data are equally spaced, that is, without loss of generality Wi,hi = 1, we obtain 



/-i 
i 



i 



\ 



D 



(8) 



V 



Tepper and Sapiro/Ll Splines for Robust, Simple, and Fast Smoothing of Grid Data 



3 



An eigendecomposition of D yields D = UAU T , where A is a diagonal matrix containing the eigenvalues of 
D, given by [20] 

Ai f-2 + 2cos((i-l)7r/n), ifi = j; 
l ' J [0, otherwise. 

Since U is a unitary matrix, we can write Equation (7) as 

y = U(I + sA 2 y 1 U T y. (10) 

Let us define the matrix r = (I + sA 2 )~ x . Trivially, 

r ■ = f[l + s(-2 + 2cos((i-l)7r/n)) 2 ]" 1 , i£i = j; 
l " 3 1 0, otherwise. 

Following Strang [15] and Garcia [7], let us observe that U T is a DCT-II matrix and U is an inverse DCT-II 
matrix. Then, Equation (10) can be expressed as 

^ = DCT- 1 (rDCT(^)), (12) 

where DCT(-) and DCT _1 () stand for the DCT-II and inverse DCT-II functions. Equation (12) provides a 
fast and simple algorithm for computing L2 splines. 

2.1. Robust estimation. 

Often in practice there are in y some values yi that could not be observed (or recorded) for some reason. 
We would like to be able to handle such cases in such a way that the missing values are inferred from the 
ones that can be observed. Let W be an n x n diagonal matrix such that W^i represents a weight assigned 
to observation i. W is defined by 

JO if datapoint i is missing; 
I p otherwise. 

where p is some arbitrary constant in (0, 1]; in practice, and without loss of generality, we set p = 1. We can 
then solve 



y — argmm 



W 1 / 2 (z-y)f 2 + s\\Dz\\ 2 2 , (14) 



which will simply omit the missing points from the computation of the residual while the regularizer will still 
have a smoothing effect over both present and missing points. Equation (14) acts as an impainting algorithm, 
filling the missing values in such a way that continuity between filled values and smoothed ones is preserved. 
The minimization of Equation (14) gives 

(I + sD T D)y = W(y-y) + y. (15) 

This leads to the iterative procedure 

f + l = (/ + SjD T D) -l ^ ( y _ y ^ + y kj ^ (16) 

which, similarly to Equation (12), becomes 

y k+1 = DCT -1 (r DCT (W (y - y k ) + y k )) . (17) 

On a different note, real data often present observations that lie abnormally far from their "true" value, 
i.e., that do not appear to follow the pattern of the other data points. The main drawback of the penalized 
least squares formulation Equation (6) is its sensitivity to these outliers. To address this issue, weights can 
be assigned to every point, as in Equation (14), such that outliers exert less influence during the estimation 
process. In this case, the weights are iteratively refined during the estimation process using robust estimators 
for the mean and variance of the data. Defining these estimators is a complex problem by itself. For details 
about how W can be set and updated for added robustness to outliers, refer to Garcia's work [7]. 
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3. LI splines 

In this section we introduce a different splines formulation in order to handle outliers in the data. We assume 
that the noise r in Equation (1) follows a Laplace distribution with zero mean and unknown scale parameter, 
which leads to K y (z) = \\z — yW-^. The regression then becomes 

min Wz-y^ + sWDzWl. (18) 

Z 

Goldstein and Osher [8] proposed a very elegant and efficient algorithm for solving the LI constrained 
problem (related to a number of very efficient optimization algorithms, e.g., see [6]) 

min 11^(^)11-, +H(V). (19) 

u 

For this, they consider the equivalent problem 

min||d|L + H(u) s.t. d = $(u), (20) 

u 

which they first convert it into the unconstrained problem 

min||d|| 1+ H( U ) + |||d-$( U )||2. (21) 

In this form, the penalty function does not accurately enforce the constraint for small A. The constraint 
is enforced by letting A —> oo. However, another solution for this new formulation is found by using the 
following two-phase algorithm 

(u k+1 ,d k+1 ) = argmin ||rf|| x + H(u) + § \\d - *(tx) - b k \\ 2 , (22) 

u 

b k+1 =b k + {${u k+1 )-d k ). (23) 

This algorithm is often denoted in the literature as split-Bregman iteration. This class of algorithms has 
several nice theoretical properties and has successfully been applied to several problems in practice such as 
image restoration [13], image denoising [18], compressed sensing [19], and image segmentation [9]; see also [6] 
and references therein. 

We use this technique for solving Equation (18). We begin by setting <S>(z) = z — y and H(z) = s \\DzW2, 
which leads to the problem 

min \\d\li +s \\Dzf 2 s.t. d = z-y. (24) 

z,d 



We then transform it into the unconstrained form 



min \\d\\ 1+ s\\Dz\\l + l\\d-z + y\\ 2 2 , (25) 

z,d 

and the Bregman iteration simply takes the form 

(z k+1 ,d k+1 ) = argmin \\dWj_ + s \\Dzf 2 + f \\d - z + y - b k f 2 , (26) 

z,d 

b k+l =}) k + ( z k+l _ y _ d k+ly ( 27 ) 

Because of the splitting of the LI and L2 components in the functional (26), we can perform this mini- 
mization efficiently by iteratively minimizing with respect to z and d separately, 

z k+1 = argmin s \\Dzf 2 + f \\d k - z + y - b k \\l , (28) 

Z 

d k+1 = argmin \\d\\ x + | \\d - z k+1 + y - b k II l . (29) 

d 

For minimizing Equation (28) we set y = d k + y — b k and s = 2s/X. We obtain 

min s\\Dz\\ 2 2 + \\y-z\\l. (30) 
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This is a classical L2 spline and can be minimized using Equation (12), as already explained. The optimal 
value of d in Equation (29) can be explicitly computed using shrinkage operators, 



d k+i 



argmm 

d 



1 1 2 
,k+l n . i uk 



x ^d-z k+1 +y-b k 



where 



= Shrink(^ +i - y + b k ,l/\), 

/ shrink(^i, 7) \ 



(31) 



Shrink(v, 7) 



shrink (u^, 7) 



\shrink(v m ,7)/ 



(32) 



and 

x 

shrink(x,7) = — max(b| —7,0). (33) 
\x\ 

We thus obtain a very efficient algorithm for computing LI splines, combining DCT and shrinkage operators. 
On a different note let us mention that Equation (25) can also be interpreted [10] as a relaxation of 

mm\\d\\ + s\\Dzf 2 + l\\d-z + y\\l, (34) 

z,d 

where the L0 norm is replaced by its (convex) LI counterpart. In this case the underlying model for y 
is y = y + r + d, where r is zero-mean Gaussian noise and d represents the "oulier" noise. Under this 
assumptions, d practically becomes an "indicator function" of the presence of (sparse) outliers (see also [3]). 
Besides the different angle in the derivation of the model, our approach differs from [10] in two very important 
points. First, we use split-Bregman iteration by introducing the variable b k in the optimization procedure, 
see equations (26) and (27). In [10] Equation (25) is first solved using direct alternate minimization over z 
and d, and then Equation (34) is solved via non-convex minimization using the previous solution as a starting 
point. Second, considering the grid structure, we use the DCT approach to solve Equation (28), instead of 
the classical Cholesky decomposition. The combination between split-Bregman and DCT results in a sound 
and fast algorithm for computing LI splines on grid data. 



3.1. Handling missing data 

In the classical L2 formulation, a diagonal binary weighting matrix W is used to cope with missing values 
(see Section 2.1 for details). Let us denote by w the diagonal of W. Let us first define 



INIi,™ = 1**1 and 

2 = 1 
Wi = l 

We then pose Equation (14) as 

V = argmin \\(z - y)\\\ w + s \\Dz\\\ , (36) 

z 

and equivalently extend Equation (18) as 



1/2 



\2,w 



E 

\ i=i 

\Wi = l 



(35) 



y = argmin \\(z - y)\\ hw + s \\Dz\\l . 

Z 



(37) 
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This leads to 

min ||d|| + a \\Dz\\l + l\\d-z + y\\ 2 2 . (38) 

z,d ' ' 

Then the split-Bregman iteration can be written as 

z k+1 = argmin s \\Dz\\ 2 2 + f \\d k - z + y - b k \\l , (39) 

Z ' 

d k+1 = argmin \\d\\ 1}W + ±\\d- z k+1 + y - b k \\l w , (40) 

d ' 
b k+l = ft fc + _ y _ d k+ly (41) 

Equation (39) can be solved using Equation (17). Solving Equation (40) amounts to performing a shrinkage 
operation on the dimensions where w equals 1. In Equation (41), it suffices to update the dimensions of b k 
where w equals 1. 



3.2. Handling multidimensional data 

Let us now return to the general case of m-dimensional data. Following Garcia [7], we extend Equation (12) 
as 

y = DCT- 1 (r m oDCT m (y)), (42) 

where DCT m (-) and DCT~ 1 (-) stand for the m-dimensional DCT-II and inverse DCT-II functions, and o 
denotes the Schur (element-wise) product. Notice that the multidimensional DCT is simply a composition of 
one-dimensional DCTs along each dimension. Extending Equation (11), r m is an m-th order tensor defined 

by 

r m = l m + (l m + s A m o A m ) , (43) 

where l m is an m-th order tensor of ones, and -f- denotes the element-wise division. Finally, A m is an m-th 
order tensor, defined by 

^,..,^=E(-2 + 2co S ^-^). (44) 

3=1 V 3 J 

where nj denotes the size of A m along the j-th dimension. 

The algorithm and its complexity. The pseudocode for the general m-dimensional case is presented 
in Algorithm 1. Let us analyze its complexity. The DCT and inverse DCT require 0(n log n) operations, 
where n = rii<j<m n r ^ ne remaining operations are linear in m. The overall complexity of the algorithm 
is then O (N Ni(m + nlogn)), where N Q and N{ are, respectively, the number of outer- loop and inner-loop 
iterations in Algorithm 1. Notice that Goldstein and Osher [8] recommend to perform only one inner-loop 
iteration for achieving optimal efficiency. Thus, we set Ni = 1 for all experiments. We will later see that in 
many cases the algorithm converges quickly (N Q can be very small then). The algorithm's complexity is thus 
dominated by the computation of the DCT and inverse DCT. Of course, these standard operations can be 
easily computed using GPU, speeding-up the execution by several orders of magnitude. 



4. Experimental results 



For all experiments we adhere to the following setup: 

1. using generalized cross validation, we find the best estimate s for s for the robust L2 formulation 
(problem (14)); 

2. we then find L2 splines (problem (6)), robust L2 splines (problem (14)), 1 and/or LI splines (prob- 
lem (18)), setting s = s. 

1 Code available at http: //www.mathworks . com/mat labcentral/f ileexchange/25634-robust- spline- smoothing- for- 1- d- 
to-n-d-data. 
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function LlSPLiNE(y, s, A, e, Ni) 

compute r m according to equations (43) and (44). 

d 1 <- 0, b 1 <- 0, k «- 1 

repeat 



for i = 1 to Ni do 



z k+i ±_ DCT- 1 (r m o DCT m (d k + y - b k )) 
d k+1 <- Shrink(z fe+1 - y + b k , 1/A) 



// usually Ni = l. 

II o denotes the element-wise product, 
//as defined in Equation (33) 



end for 



fc <- fc + 1 



until -z k \\ 2 I \\z k \\ 2 >e 



return z' 



end function 



Fig 1. Pseudocode of the LI spline for evenly spaced-data via Bregman iteration 



This protocol allows us to show that, even when s is chosen to fit optimally the robust L2 formulation, 
the proposed method provides better estimates. For the LI formulation, in Equation (18), we simply set 
A = min(s, 1) for all examples. We recall that Ni is set to 1. We also set e — 10 -3 (see Algorithm 1) and 
additionally limit the maximum number of outer iterations to a hundred. The algorithm stops when any of 
the two conditions is met. 

Fig. 2 presents two one-dimensional examples. We depict the original signal y G [1, . . . , n] —> M, where 
n = 2 16 . We observe the signal y = y + ri, where r\ is Gaussian noise. Some points yj (1 < j < n) are further 
contaminated with uniform noise where G [a, . . . , 6], such that yj = min(max(^j + r\ + r2, a), b). The 
points affected by T2 are depicted in red and the remaining ones in green. In the top row, a = — 5, b = 5; and 
in the bottom row, a = 0, b = 5. In both examples, only the LI spline is correct. The classical and robust L2 
splines are both unable to correctly recover the original data in the corrupted part. 

Fig. 3 shows the evolution of the relative error as the number of iterations increases for the example in 
Fig. 2 (top row). As we can observe, the proposed algorithm is able to converge quickly, reaching a precision 
of 10 -3 in less than 20 iterations. 

Fig. 4 depicts the relative time-cost of each operation during the execution of the proposed method 
(Fig. 2, top row). Computing the DCT and the inverse DCT covers more than 84% of the total running 
time. Implementing these standard operations in GPU would boost the performance of the algorithm by 
orders of magnitude. 

In Fig. 5 we present a two-dimensional example. We depict the original signal y G [1,...,256] 2 — » 
[—6.5497,8.1054] in Fig. 5a, and we add two types of noise: first, Gaussian noise r\ with zero mean and 
variance a 2 = 2 (Fig. 5b), and then uniform noise T2 in the interval [—5 • max(y + r\ ),..., 5 • max(y + n)] 
(Fig. 5c). Again, only the LI spline correctly recovers the original signal. 

We next test the proposed algorithm with a climate time-series provided by the Met Office Hadley Cen- 
tre [4]. 2 The dataset contains the evolution of global average land temperature anomaly (in °C) with respect 
to the 1961-1990 average temperature. The results, which confirm an upward trend in the second half of the 
20th century, are shown in Fig. 6. 

We also test on this dataset the effect of varying the parameters s and A, see Fig. 7. As in the classical 
L2 formulation, s has a direct impact on the obtained result, controlling the degree of smoothness of the 
solution, see Fig. 7a. On the contrary, Fig. 7b shows that the newly introduced parameter A is very stable 
and provides very similar results in a wide range (A G [0.1, 100]). This stability allows us to fix its value to 
A = 1 for all the experiments in this work. 

Another interesting example is presented in [10]. The dataset consists of power consumption measurements 
(in kW) for a government building, collected every fifteen minutes from July 2005 to October 2010. As in [10], 
we downsample the data by a factor of four, yielding one measurement per hour, and use only a subset of 
the whole data. The results are displayed in Fig. 8. 

We also test in a synthetic example the ability to recover signals with sharp transitions, see Fig. 9. In 
this case we use a simple piece-wise constant function. We can observe clear overshoot (plus ringing) effects 
on the L2 and robust L2 splines. The robust L2 spline also results in transitions with less vertical slopes, 



2 Data are available in http: //hadobs .metof f ice . com/crutem3/diagnostics/global/nh+sh/annual. 
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(a) L2 spline (b) Robust L2 spline (c) LI spline 



Fig 2. Noisy data (in green) further contaminated with uniform noise (in red). The spline y, depicted in blue, is obtained using 
different fitting terms. The proposed method with the LI fitting term recovers the correct shape. 



0\ 




Iterations 

Fig 3. The relative error (in logarithmic scale) as the iterations progress when computing the first example in Figure 2. The 
algorithm is able to quickly decrease the error during the first 20 iterations. 
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Fig 4. Percentage of the execution time spent in each operation when computing the first example in Figure 2. Clearly, the 
vast majority of time is spent in DCT or inverse DCT operations. 




(d) L2 spline (e) Robust L2 spline (f) LI spline 



Fig 5. Synthetic two-dimensional example: the original data y is contaminated with Gaussian noise r\ and then with a large 
uniform noise r<i. With this input signal, y + r\ + r<i, the proposed method is the only one able to recover the correct shape. 
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lr 

Original signal 




1840 1860 1880 1900 1920 1940 1960 1980 2000 2020 

Years 

Fig 6. Global average land temperature anomaly (°C) with respect to 1961-1990: smoothed versus original year-averaged data [4]. 
Control points are points such that yi = yi . 



creating a bluring effect. With the LI spline we obtain a much better reconstruction, with almost non-existent 
overshooting. 

This very same effect can be observed in real examples, see Fig. 10. When approximating images with 
splines, some structure is lost by blur and some structure is artificially created by overshooting and ringing. 
This can be observed in Figs. 10a and 10b, were the difference between the original image and the image 
estimated by robust L2 splines exhibits structure. Observe, however, that almost no structure in the difference 
is visible when the reconstruction is performed using LI splines. 

J^.l. Application to range data 

In this section we perform smoothing of depth data obtained with a Kinect camera. This kind of data is 
particularly challenging because: 

• it presents relatively smooth areas separated by sharp transitions, 

• edges are highly noisy, that is, edge pixels oscillate over time between foreground and background, and 

• it contains missing data, which appear for two different reasons: (1) the disparity between the IR 
projector and the IR camera produces "shadows," and (2) the depth cannot be recovered in areas 
where the IR pattern is not clearly observable (e.g., because they receive direct sunlight or interference 
from another Kinect). 

We use splines to interpolate and denoise these data, showing the advantage of LI splines over its robust L2 
counterpart. The displayed images are part of the LIRIS human activities dataset [17]. 

In the first example, shown in Fig. 11, we use a single depth frame (with standard Kinect resolution 
of 640 x 480). The missing data are represented in black, while depth data goes from red to yellow as 
depth increases. Both, the LI spline and the robust L2 spline are able to interpolate the missing data with 
reasonable values. Notice, however, that the latter exhibits, as aforementioned, overshooting and ringing 
(clearly perceived in the ID profile). These effects are much milder in the LI reconstruction. 

In Fig. 12a we can clearly observe that the position of the missing data is not consistent across frames. 
We can integrate data from several frames to achieve more accurate interpolations, by performing 3D re- 
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Original signal 

■ s=100 

■ s=200 

■ s=300 
s=400 




1840 1860 1880 



1900 

(a) Variable s, fixed A = 1 



1920 1940 1960 
Years 



1980 2000 2020 




1840 1860 1880 1900 1920 1940 1960 1980 2000 2020 

Years 



(b) Fixed s = 200, variable A. 



Fig 7. Changing the value of s in LI splines controls the trade-off between fitting and smoothing. LI splines are quite insensitive 
to the value of \, making it easy to tune. 
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Fig 8. Power consumption measurements (in kW) for a government building [10] with zoom-in details on the bottom. 
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(a) L2 spline 
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(b) Robust L2 spline 



(c) LI spline 



Fig 9. In the case of a periodic piece-wise constant signal, the L2 and robust L2 approximations exhibit overshoot and ringing 
(see zoom-in details on the bottom). These effects are much attenuated by the LI spline approximation. 



Original image 



Noisy image 





Robust L2 spline 



LI spline 
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(a) 



(b) 



Fig 10. (a) When no noise is added, the robust L2 spline removes much more structure than the LI spline. For the robust L2 
spline and for the LI spline, the norm of the difference is 9.9866 and 0.7301, respectively, (b) When salt & pepper noise is 
added, the robust L2 spline removes noise and structure, while the LI spline removes much less structure but leave some noise. 
Notice that we are not proposing a new image- denoising algorithm, we just use images for showing that LI splines respect more 
the signal structure than L2 splines. 
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(a) Original depth image 



(b) Robust L2 spline 



(c) LI spline 





Original values 
- Robust L2 spline 
-LI spline 



200 300 400 500 600 700 100 200 300 400 500 

(d) Profile of a single row. Left, original values; right, original and reconstructed values. 



Fig 11. Smoothing and interpolating using a single depth frame. We simply fit a 2D spline to the depth image. Since the data 
presents several "jumps/ 7 the robust L2 spline must under-smooth the data to be able to fit it correctly. The LI spline presents 
a good trade-off between fitting and smoothing. 
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Table 1 

Execution times (in seconds) and number of iterations until convergence of the proposed algorithm for the different 

experiments performed in this work. 









Size 


Robust L2 spline 


LI spline 








Time 


Time 


Iters. 




Fig. 


2 


2 20 


4.931 


3.590 


7 




Fig. 


5 


256 x 256 


0.541 


0.982 


100 


Complete 
data 


Fig. 
Fig. 


6 

8 


163 
501 


0.007 
0.010 


0.045 
0.010 


72 
11 


Fig. 


9 


2 io 


0.035 


0.007 


6 




Fig. 


10a 


321 x 481 


1.167 


0.758 


17 




Fig. 


10b 


321 x 481 


1.143 


0.873 


20 



Missing Fig. 11 480 x 640 0.911 0.266 2 

data Fig. 12 480 x 640 x 3 7.511 2.202 3 



constructions. Thus, in this example, we treat depth data as a 3D signal (2D + time), by considering three 
consecutive frames. The data dimensionality is then 640 x 480 x 3. A full depth video can be smoothed by 
using 3D splines as a sliding-window type of filter. The robust L2 spline again presents a noisier behavior 
and with significative overshooting. On the other hand, the LI spline is much smoother in smooth areas 
while correctly preserving abrupt transitions. 

Running times. We present in Table 1 the running-time and number of iterations until convergence for 
every example in this work. The time of the robust L2 and the LI splines is comparable. All code is written 
in pure Matlab, with no C++ or mex optimizations. All experiments were run on a MacBook Pro with a 
2.7GHz Intel Core i7 processor. Finally, note that in most cases the algorithm converges in less than twenty 
outer- iterations (recall that e = 10 -3 ). In the example in Fig. 5, the maximum number of iterations (100) is 
reached with a final error of 10 -2,8 . 

5. Conclusions 

We have presented a new method for robustly smoothing regularly sampled data. We do this with modified 
splines, where we replace the classical L2-norm in the fitting term by an Ll-norm. This automatically handles 
outliers, thus obtaining a robust approximation. 

We also presented a new technique, using split-Bregman iteration, for solving the resulting optimization 
problem. The algorithm is extremely simple and easy to code. The method converges very quickly and has 
a small memory footprint. It also makes extensive use of the DCT, thus being straightforward to implement 
in GPU. These characteristics make this method very suitable for large-scale problems. 
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